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Abstract — We consider finite-dimensional Markovian open 
quantum systems, and characterize tlie extent to wliich time- 
independent Hamiltonian control may allow to stabilize a target 
quantum state or subspace and optimize the resulting con- 
vergence speed. For a generic Lindblad master equation, we 
introduce a dissipation-induced decomposition of the associated 
HUbert space, and show how it serves both as a tool to 
analyze global stability properties for given control resources 
and as the starting point to synthesize controls that ensure 
rapid convergence. The resulting design principles are illustrated 
in realistic Markovian control settings motivated by quantum 
information processing, including quantum-optical systems and 
nitrogen-vacancy centers in diamond. 



I. Introduction 

Devising effective strategies for stabilizing a desired quan- 
tum state or subsystem under general dissipative dynamics 
is an important problem from both a control-theoretic and 
quantum engineering standpoint. Significant effort has been 
recently devoted, in particular, to the paradigmatic class of 
Markovian open quantum systems, whose (continuous-time) 
evolution is described by a quantum dynamical semigroup [1 1. 
Building on earlier controllability studies ||2|, jS), Markovian 
stabilization problems have been addressed in settings ranging 
from the preparation of complex quantum states in multipartite 
systems to the synthesis of noiseless quantum information 
encodings by means of open-loop Hamiltonian control and 
reservoir engineering as well as quantum feedback Q, jS), jBj, 
El, m, E). While a number of rigorous results and control 
protocols have emerged, the continuous progress witnessed by 
laboratory quantum technologies makes it imperative to de- 
velop theoretical approaches which attempt to address practical 
constraints and limitations to the extent possible. 

In this work, we focus on the open-loop stability properties 
of quantum semigroup dynamics that is solely controlled in 

Francesco Ticozzi is with the Dipartimento di Ingegneria dell'Informazione, 
Universita di Padova, via Gradenigo 6/B, 35131 Padova, Italy (email: 
ticozzi@dei.unipd.it). 

Riccardo Lucchese is with the Dipartimento di Ingegneria 
dell'Informazione, Universita di Padova, via Gradenigo 6/B, 35131 
Padova, Italy (email: lucchese@dei.unipd.it). 

Paola Cappellaro is with Department of Nuclear Science and Engineering, 
Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, 
MA 02139, USA (email: pcappell@mit.edu). 

Lorenza Viola is with Department of Physics and Astronomy, Dart- 
mouth College, 6127 Wilder Laboratory, Hanover, NH 03755, USA (email: 
lorenza.viola@dartmouth.edu). 

F. T. acknowledges hospitality from the Physics and Astronomy Department 
at Dartmouth College, where this work was performed, and support from 
the University of Padova under the QUINTET project of the Department of 
Information Engineering, and the QFUTURE and CPDA080209/08 grants. 



terms of time-independent Hamiltonians, with a twofold mo- 
tivation in mind: (i) determining under which conditions a de- 
sired target state or subspace may be stabilizable given limited 
control resources; (ii) characterizing how Hamiltonian control 
influences the asymptotic speed of convergence to the target 
set. A number of analysis tools are obtained to this end. On 
the one hand, we introduce the notion of a dissipation-induced 
decomposition of the system Hilbert space as a canonical 
way to represent and study Markovian dissipative dynamics 



(Sec. III-A I. A constructive algorithm for determining such a 



(unique) decomposition is presented, as well as an enhanced 
algorithm that finds which control inputs, if any, can ensure 
convergence under fairly general constraints (Sec. |I1I-B[ ). On 
the other hand, we present two approaches to analyze the 
convergence of the semigroup to the target set: the first, which 
is system-theoretic in nature, offers in principle a quantitative 
way to computing the asymptotic speed of convergence (or 
Lyapunov exponent) from a certain reduced dynamical matrix 
(Sec. IV-A the second, which builds directly on the above 



dissipation-induced decomposition and we term connected 
basins approach, offers a qualitative way to estimating the 
convergence speed and designing control in situations where 
exact analytical or numerical methods may be impractical 
(Sec. |1V-B| |. By using these tools, we show how a number 
of fundamental issues related to role of the Hamiltonian in 
the convergence of quantum dynamical semigroups can be 
tackled, also leading to further physical insight on the interplay 
between coherent control and dissipation |10|. A number of 
physically motivated examples are discussed at length in Sec. 
[V] demonstrating how our approach can be useful in realistic 
quantum control scenarios. 

II. Quantum dynamical semigroups 

A. Open-loop controlled QDS dynamics 

Throughout this work, we shall consider a finite- 
dimensional open quantum system with associated complex 
Hilbert space H, with dim('H) — n. Using Dirac's notation 
[11|, we denote vectors in "H by and linear functionals in 
the dual "H^ ~ "H by {4>\. Let in addition S('H) be the set of Hn- 
ear operators on T-L, with t}{H) being the Hermitian ones, and 
2) (7^) C 1)(H) the trace-one, positive semi definite operators 
(or density operators), which physically represent the states of 
the quantum system. The class of dynamics we consider are 
one-parameter quantum dynamical semigroups (QDSs), that 
is, continuous families of trace-preserving completely positive 
maps {7t}t>o that enjoy a forward (Markov) composition law. 
The evolution of states p G ©("H) is then governed by a master 
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equation Ea, IH, IHI, Q of the Lindblad form: 

p , 
P = C[p] = -j \H,p\+Y, (l.pLI - -{LlLk,p}) , 



(1) 



fe=i 



where the effective Hamihonian H £ ^nd the noise 

operators {Lfe} C B{H) describe, respectively, the coherent 
(unitary) and dissipative (non-unitary) contributions to the 
dynamics. In what follows, h ^ 1 unless otherwise stated. 

We focus here on control scenarios where the Hamiltonian 
H of the system can be tuned through suitable control inputs, 
that is. 



Hiu) = Ho 



i=i 



(2) 



where Hq = Hq represents the free (internal) system Hamilto- 
nian, and the controls u, G M modify the dynamics through the 
Hamiltonians Hj = . In particular, we are interested in the 
case of constant controls Uj taking values in some (possibly 
open) interval Cj C M. The set of admissible control choices 
is then a subset C CR". 

B. Stable subspaces for QDS dynamics 

We begin by recalling some relevant definitions and results 
of the linear-algebraic approach to stabilization of QDS de- 
veloped in m, 0, IS), nS], im. consider an orthogonal 
decomposition of the Hilbert space H := Hs ® Hb., with 
dim('H5') — m < n. Let {|si)} and be orthonormal 

sets spanning T-Ls and T-Lb, respectively. The (ordered) ba- 
sis {|si), . . . , |s„i), |ri), . . . , |r„_„i)} induces the following 
block structure on the matrix representation of an arbitrary 
operator X G Bin): 

X=\f (3) 

Let in addition supp(X) := kcr(X)-'-. It will be useful to 
introduce a compact notation for sets of states with support 
contained in a given subspace: 

'ps 0" 




35(H) :={pgS)(H) 



P 



, ps€^{Hs)}. 



As usual in the study of dynamical systems, we say that 
a set of states W is invariant for the dynamics generated by 
C if arbitrary trajectories originating in >V at ^ = remain 
confined to W at all positive times. Henceforth, with a slight 
abuse of terminology, we will say that a subspace Hs C H 
is C-invariant (or simply invariant) when 3s{'H) is invariant 
for the dynamics generated by C. 

An algebraic characterization of "subspace-invariant" QDS 
generators is provided by the following Proposition (the proof 
is given in ID, in the more general subsystem case): 

Proposition 1 ( S -Subspace Invariance): Consider a QDS 
on "H = 'Hs®'Hr, and let the generator £ in ([TJ be associated 
to an Hamiltonian H and a set of noise operators {Lk}. Then 
is invariant if and only if the following conditions hold: 



Lp.k = 0, 



k 

Ls.k 




P,k 
R,k 



(4) 



The following Corollary can be derived in a straightforward 
way from Proposition [T] and it will be key to establishing the 
results of the next section: 

Corollary 1 (R-Subspace Invariance): Consider a QDS on 
H — Hs ® Hji, and let the generator £ in ([T]l be associated 
to an Hamiltonian H and a set of noise operators {Lk}- Then 
Hr is invariant if and only if the following conditions hold: 



fe 

'Ls,k 
Lg^k Lp^k 



Lk — 



(5) 



V/c. 



One of our aims in this paper is to determine a choice con- 
trols that render an invariant subspace also Globally Asymp- 
totically Stable (GAS). That is, we wish the target subspace 
Hs to be both invariant and attractive, so that the following 
property is obeyed: 



lim S{Tt{p),3s{H)) = 0, 

t^OO 



VpeD(H), 



where S{a, W) := inirew ll"'^'''!!- In lU, a number of results 
concerning the stabilization of pure states and subspaces by 
both open-loop and feedback protocols have been established. 
For time-independent Hamiltonian control, in particular, the 
following condition may be derived from (|4]) above: 

Corollary 2 (Open-loop Invariant Subspace): Let H = 
'Hs®T~Lr. Assume that we can modify the system Hamiltonian 
as H' = H+Hc, with He being an arbitrary, time-independent 
control Hamiltonian. Then Jsi'H) can be made invariant under 
£ if and only if Lq k = for every k. 

The following theorems (proved in 15]) provide a general 
necessary and sufficient condition for attractivity, and charac- 
terize the ability of Hamiltonian control at ensuring the desired 
stabilization: 

Theorem 1 (Subspace Attractivity): Let % = "Hs ® 'Hr, 
and assume that Hs is an invariant subspace for the QDS 
dynamics in Q. Let 



Pi ker(L 



P.k) 



(6) 



k=l 



with each matrix block Lp k representing a linear operator 
from Hr to Hs- Then Hs is GAS under C if and only if 
Hr' does not support any invariant subsystem. 

Theore}7i 2 (Open-loop Subspace Attractivity): Let H = 
Hs © Hr, with Hs supporting an invariant subsystem. As- 
sume that we can apply arbitrary, time-independent control 
Hamiltonians. Then 3s{'H) can be made GAS under £ if and 
only if 3r{H) is not invariant. 

From a practical standpoint, a main limitation to imple- 
menting the constructive procedure developed in the proof of 
Theorem |2] is that the assumption of access to an arbitrary 
control Hamiltonian He is typically too strong. Thus, we shall 



develop (Sec. III-B 1 an approach that allows us to characterize 
whether and how a given stabilization task can be attained with 
available (in general restricted) time-independent Hamiltonian 
controls, as well as the role of the Hamiltonian component in 
determining the speed of convergence. 
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III. Analysis and synthesis tools 

A. Dissipation-induced decomposition and stability 

Suppose that we are given a target subspace Hs C H. By 
using Propositions [T] one can easily check if it is invariant 
for a given QDS. If this is the case, we next exploit the ideas 
of Theorem [T] to devise a constructive procedure that further 
determines whether Hs is attractive or not. Explicitly, this can 
be attained by constructing a special (canonical) Hilbert space 
decomposition into orthogonal subspaces, which is induced by 
the matrices associated to H, {L^} in ([TJ. 

Algorithm for Dissipation-Induced Decomposition 

Assume to be invariant. Call "H^"* := T-Lr, := Us, 
choose an orthonormal basis for the subspaces and write the 
matrices with respect to that basis. Rename the matrix blocks 



(0) 



(0) 
'R.k 



(0) 
^S.k ■~ 



For j > 0, consider the following iterative procedure: 

1) Compute the matrix blocks Lp'f, according to the de 



composition H'-^'^ ^ Ug ® U 

2) Define flfckerL^J^. 

3) Consider the following three sub-cases: 



a. IfH 



(j+i) _ 



R 



= {0}, define := U 



R ■ 



The iterative procedure is successfully completed. 



b. Ifn^d^'^ C define n^^+'^ := H^^eH^+'^ 



c. If = n'^' (that is, L'^'^ = Vfc), define 



(3) 



2 ^Q,k^R,k- 
k 



- If ^ 0, re-define H^"^^^ := kcr(£^''). 



H^' and 



If 7^(^+1) = {0}, define % 
the iterative procedure is successfully completed. 
Otherwise define H^^^' := "H^^ v!^^^^^- 



- If £ 



0, then, by Corollary ^ 'H}^' is 



invariant, and thus (by Theorem [TJ "Hs cannot 
be GAS. 



at i = q. Then we have obtained a decomposition 1-Lr = 



(2). 



and we can prove by (finite) induction 



that no invariant subspace is contained in Hr. 

Let us start from 'h!j^\ By definition, since the algo- 



rithm is completed when H 



nfcker(L^i) 

Theorem 111 and Corollary 111 not only does Wj!' fails to be 
invariant, but it cannot contain an invariant subspace. 

Now assume (inductive hypothesis) that (S- ■ ■ 

£ + 1 < q, does not contain invariant subspaces, and that 



{0} or 



(i+i) 

R 



n 



(9+1) 
R 



{0}, either 



is full-rank. In either case, by 

(<?) 



(by contradiction) K 



does. Then 



the invariant subspace should be non-orthogonal to Hj^ , 
which is, by definition, the orthogonal complement of either 



Pl/j ker(Lp^.^'') or ker(£p '''). But then any state p with 
support only on H 



n 



(^+1) 



n^r^^ and non-trivial 



support on Wrp would violate the invariance conditions and 

(i) ' ' 

"exit" the subspace. Therefore, Wrp ® . . 



"HI?^ does not 



contain invariant subspaces. By iterating until I' = 1, we have 
that T-Lr cannot contain invariant subspaces and, by Theorem 
[T] the conclusion follows. ■ 
Formally, the above construction motivates the following: 
Definition 1: Let '3s{T~L) be GAS for the dynamics ([T]i. The 
Hilbert space decomposition given by 



(1) 



(2) 



(«) 



(7) 



as obtained with the previous algorithm, is called Dissipation- 
Induced Decomposition (DID). Each of the subspaces H^j) in 
the above direct sum (or possibly a union of them) will be 
referred to as a basin. 

Partitioning each matrix in blocks according to the 
DID results in the following canonical structure, where the 
upper block-diagonal blocks establish the dissipation-induced 



connections between the different basins Ti 





' Ls 


f (0) 

i^p 








Lrp 


L^^^ ... 


Lk — 






r(2) f(l) .. 

Ijrp ^ P 



(8) 



4) Define U 



(j+i) a/0) 

— 7X0 



f 



(i+i) 



To construct a basis The blocks D^'fj, are the restrictions of i^. to the subspaces 



for T-Lg^'^' , append to the already defined basis for ^-t' ■ Since in step 3.b of the DID algorithm H^rji' is defined 



,U) 



(i+1) 



an orthonormal basis for "H 
5) Increment the counter j and go back to step 1). 



The algorithm ends in a finite number of steps, since at every 

(i) 

iteration it either stops or the dimension ofHp is reduced by 
at least one. As anticipated, its main use is as a constructive 
procedure to test attractivity of a given subspace T-Ls- 

Proposition 2: The algorithm is successfully completed if 
and only if the target subspace Hs is GAS (JsiH) is GAS). 

Proof: If the algorithm stops because Cp = for some 
j, then Corollary [T] implies that T-Lr contains an invariant 
subspace, hence Hs cannot be GAS. On the other hand, let 
us assume that the algorithm runs to completion, achieved 



to be in the complement of Hp^^^ = ker Lp],., at each 
iteration the only non-zero parts of the Lp blocks must be in 
the + 1) block, which we have denoted by L^p\ in ([Sj. 
In the upper-triangular part of the matrix, the other blocks 
of any row are thus zero by construction. If some Lp\ — 
V k, then either the dynamical connection is established by the 
Hamiltonian H, through the block Hp^ (as checked in step 
3.c), or the target subspace is not GAS. 

Remark: The DID in (|7| is unique, up to an arbitrary choice 
of basis in each of the orthogonal components Hs, H-i^, i — 
1,. . . ,q. This freedom may be exploited to further put the 
matrix blocks in some (problem-dependent) convenient form. 

We illustrate the DID algorithm with an explicit Example, 
which will also be further considered in Sec. |V] 
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Example 1: Consider a bipartite quantum system consisting 
of two two-level systems (qubits), and on each subsystem 
choose a basis {|0)„, |l)ri}, with n = 1,2 labeling the qubit. 
The standard (computational) basis for the whole system is 
then given by {|00), |01), |10), |11)}, where \xy) := \x)i <g> 
\y)2- As customary, let in addition {ctq, a = x,y,z} denote 
Pauli pseudo-spin matrices ifTTIl . with the "ladder" operator 
a_|_ := {ax + i(Jy)/2 = |0)(1|. Assume that the dynamics is 
driven by the following QDS generator: 



= C[p] = -t[H, p] + LpL^ - ^{L^L, p}, 



(9) 



where 



H = il(T,+crx)(Sl + I®{-l(7,+ax), (10) 

L = (T+®I + (11) 



It is easy to verify that the (entangled) state p^ = |V'o)(V'oli 
with 

|^o) = 4.(|00>-|01> + |10)). 



V3 

is invariant, that is, C[pd\ = 0. We can then construct the DID 
and verify that such state is also GAS. 

By definition, H^g^ = span{|'i/'o)}, and one can write its 
orthogonal complement as — span{|V'i), |'02), IV's)}; 

with an orthonormal basis being given for instance by: 

1 



^(|01> + |10)), 
i(|01)-|10)) 



IV-l) = 111), IV^2> ^ 

l^3) = -^(|00)- 
We begin the iteration with j 

[0 0.816 0]. We move on (step 2), by defining H^^^ := 
ker(Lp^) = spanjlV'i), IV's)}- We next get (step 3.b): 



(step 1), having Lp 



(12) 

(13) 

(0) _ 



(1) 



so that (step 4): 



n 



(1) 



n 



©H^^^ =span{|V'o),|V'2)}. 



(0) 



/(I) 



span{|?/>2)}, 



We thus set j = I and iterate, obtaining: 





1.414 



n 



(2) _ 



R =kcr(i'^ 

/(2 



span{|V'3)}, y-T^ ^spanjlV-i)} 



H);> ^span{\^o)A^2)Ai^l)}■ 



ln the third iteration, with j = 2, we find that Ly = [0 0]^ 

(3) (2) 

Hence we would have Hy = 'Hjj , so we move to step 3.c 
By computing the required matrix blocks we get: 



(2) _ 



-iH 



(2) 



Re-defining H^^^ 



n 



(3) 



^ EOS = -[0-1-732 op 



ker(£p), we find that "hL^^ = {0}, thus 



(3) 

y.)^ and the algorithm is successfully completed. 



Having support on Tis alone, pd is thus GAS, and in the 
ordered basis {[^Ao), |'02), IV'i): IV's)}, consistent with the DID 



Hs ® H. 



(1) 



n 



(2) 



(3) 

Hj^ , we have the following matrix 



representations (cf. Eq. m 



L = 



■ 


0.816 





" 








1.414 




















-1.155 











' 



















1.414 


-1.732 


H = 





1.414 













-1.732 









It is thus evident how the transitions from "hI^^' to "Hs, and 



f 2 

from n)f' to n'j:', are enacted by the dissipative part of the 

f3) 

generator, whereas only the Hamiltonian is connecting Hj,' 
to 'H^\ destabiUzing \ip3)- 



/(I) 



B. QDS stabilization under control constraints 

The DID algorithm can be used as a design tool to determine 
whether an admissible control Hamiltonian may achieve sta- 
biUzation under constrained capabilities, . . . , u^) e C C 
My . Assume we are given a target H5, which need not be 
invariant or attractive. We can proceed in two steps. 

1 ) Imposing invariance: Partition H, according to H = 
Hs © Hr- If Lq^k 7^ for some k, then Hs is not invariant 
and it cannot be made so by Hamiltonian control, hence it 
cannot be GAS. On the other hand, if Lp fe — for all k, then 
Hs cannot be made GAS by Hamiltonian open-loop control 
since 3]i{H) would necessarily be invariant too (Theorem [2|l. 
When Lq,/c = for all k and there exists a k such that 
7^ 0, we need to compute (Proposition |T|i 



Cp{u) = iHp{u) 



S,k^PM- 



If Cp{u) 7^ for all u £ C, then the desired subspace 
cannot be stabilized. Let C'"^ be the set of controls (if any) 
such that if u e C^°\ then Cp{u) = 0. 

2j Exploring the control set for global stabilization: Hav- 
ing identified a set of control choices that make Hs invariant, 
we can then use the DID algorithm to check whether they can 
also enforce the target subspace to be GAS. 

By inspection of the algorithm, the only step in which a 
different choice of Hamiltonian may have a role in determining 
the attractivity is 3.c. Assume that we fixed a candidate control 
input u, we are at iteration j and we stop at 3.c. Assume, 
in addition, that the last constrained set of controls we have 
defined is C(^\ < ^ < j (in case the algorithm has not 
stopped yet, this is C^^-'). Two possibilities arise: 

. If C^p ^ 0, define C^^^ as the subset of C^^^ such that if 
u G C'^^\ then it is still true that Cp{u) ^ 0. Pick a choice 
of u e C^-'^ , and proceed with the algorithm. Notice that if 
there exists a control choice ii such that Cp{u) has full 
rank, we can pick that and stop the algorithm, having 
attained the desired stabilization. 

• If Cp = 0, the algorithm would stop since there would 



f 7 + 1 1 17 1 

be no dynamical link from H^ ' towards H j. , neither 
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enacted by the noise operator nor by the Hamiltonian. 
Hence, we can modify the algorithm as follows. Let us 
define C'^^ as the subset of C^'^^ such that if fi e C^^\ then 
Cp{u) 7^ 0. If C^^-' is empty, no other choice of control 
could destabilize H^^^^ and Hs cannot be rendered 
GAS. Otherwise, redefine C^^^ := C'^\ pick a choice 
of controls in the new C^^^ (for instance at random), and 
proceed with the algorithm going back to step I. 
The above procedure either stops with a successful com- 
pletion of the algorithm or with an empty C*^^-*. In the first 
case the stabilization task has been attained, in the second it 
has not, and no admissible control can avoid the existence of 
invariant states in TIr. 

Note that if each Cj (thus C) is finite, for instance in the 
presence of quantized control parameters, the algorithm will 
clearly stop in a finite number of steps. More generally, in the 
following Proposition we prove that in the common case of 
a multi-interval as the set of admissible controls, the design 
algorithm works with probability one: 

Proposition 3: If C is a bounded multi-interval of M.", the 
above modified algorithm will end in a finite number of steps 
with probability one. 

Proof: The critical point in attaining GAS is finding a 
set of control values that ensures invariance of the desired set 
when the the free dynamics would not. In fact, to this end we 
need to find a w e C(°) = {u G C\Cp{u) ^ 0}. Being C(°) the 
intersection between a multi-interval and a (j/— l)-dimensional 
hyperplane in W , the set C'*^-* belongs to a lower-dimensional 
manifold than C. 

Once invariance has been guaranteed, we are left with the 
opposite problem: at each iteration j, we need to avoid the 
set of controls such that £p' = 0. This is again a. [v — 1)- 
dimensional hyperplane in M". Therefore, if a certain uq is 
such that Cp = but not all of them are, this belongs to 
a lower-dimensional manifold with respect to C'"'. Hence, 
picking a random choice of w e C'"-* (with respect to a uniform 
distribution) will almost always guarantee that the algorithm 
will stop in a finite number of steps. ■ 

C. Approximate state stabilization 

A necessary and sufficient condition for a state (not nec- 
essarily pure) to be GAS is that it is the unique stationary 
state for the dynamics ID: this fact can be exploited, under 
appropriate assumptions, to approximately stabilize a desired 
pure state pd when exact stabilization cannot be achieved. 

Assume that at the first step in the previous procedure we 
see that pd is not invariant, even if ^ ~ for all fc, since 

r(o) _ 1 \^ r(0)tr(0) / n 

L,p - lUp - 2 2^ S,k ^P,k r 

k 

and there exists no choice of controls that achieve stabilization. 
If however the (operator) norm of £p' can be made small, in 
a suitable sense, we can still hope that a GAS state close to 
Pd exists. This can checked as it follows: 



Define Hp 



iC^p\ Consider a new Hamiltonian 



AH = 



'Hs 


Hp 




■ 


Hp 


Hp 


Hp 


+ 


[Hi 






By construction, pd is invariant under H. 

• Proceed with the algorithm described in the previous 
subsection in order to stabilize pd with H instead of H. 

• As a by-product, the subset of control values that achieve 
stabilization is found. Let it be denoted by 5 C C; 

• Determine e 5 such that min„g5 0p\u)\\^ IS 
attained. 

After the control synthesis, the generator for the actual 
system is in the form: 

p = C[p] ^ AC[p], 

with AC[p] = —i[AH,p], and C[p] having pd as its unique 
stationary state. It can be easily proved (by using the standard 
coherence-vector representation, following the same argument 
of Section III.E in |4|) that the perturbed generator will still 
have a unique stationary state provided that the norm of AC is 
small, with respect to the smallest absolute value of the non- 
zero eigenvalues of C. In our setting, ||A£j| can be bounded 
by twice the (operator) norm of Hp. By continuity arguments, 
one can also show that the new stationary state will then 
be close to the target one. This ensures stabilization of a 
(generally) mixed state in a neighborhood of pd or, in the 
control-theoretic jargon, "practical stabilization" of the target 
state, the size of the neighborhood depending on ||A£||. 

IV. Speed of convergence of a QDS 

How quickly can the system reach the GAS subspace Hs 
from a generic initial state? We address this question in two 
different ways. The first approach relies on explicitly comput- 
ing the asymptotic speed of convergence by considering the 
spectrum of C (sp(£) henceforth) as a linear superoperator 
Despite its simplicity and natural system-theoretic interpreta- 
tion, the resulting (worst-case) bound provides little (if any) 
physical intuition on what effect individual control parameters 
have on the overall dynamics. In the second approach, which 
builds directly on the DID, we argue that convergence is 
limited by the slowest speed of transfer from a basin subspace 
to the preceding one in the chain. Despite its qualitative nature, 
this has the advantage of offering a more transparent physical 
picture and eventually some useful criteria for the design of 
rapidly convergent dynamics. 

A. System-theoretic approach 

The basic step is to employ a vectorized form of the QDS 
generator £ (also known as superoperator or Liouville space 
formalism in the literature |161), in such a way that standard 
results on linear time-invariant (LTI) state-space models may 
be invoked. Recall that the vectorization of a n x m matrix 
M, denoted by vec(M), is obtained by stacking vertically the 
m columns of M, resulting in a -dimensional vector [TTJ. 
Vectorization is a powerful tool when used to express matrix 
multiplications as linear transformations acting on vectors. The 
key property we will need is the following: For any matrices 
X, Y and Z such that their composition XYZ is well defined, 
it holds that: 



vec{XYZ) = {Z^ (g> X)vec(y), 



(14) 
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where the symbol (g) is to be understood here as the Kro- 
necker product of matrices. The following Theorem provides a 
necessary and sufficient condition for GAS subspaces directly 
in terms of spectral properties of the (vectorized) generator 
(compare with Theorem [TJ: 

Theorem 3 (Subspace Attractivity): Let H = Us ® 'Hr, 
and assume that Hs is an invariant subspace for the QDS 
dynamics in ([T]). Then Hs is GAS if and only if the linear 
operator defined by the equation 



k 

- ijY.'^R® {L^p,kLp,k + L^R,kLR,k) (15) 

k 

- 9 'Y{L'p,kL*p.k + LR^kL*R,k) ® Ifl- 



does not have a zero eigenvalue. 

Proof: Let 11^^ = [O 1 r\ ■ explicitly computing the 
generator's i?-block we find: 

TiRClppj, - - '-{[Hr, pr] + hIpp - p^pHp) 

+ i^Q kPs + LR,kPp)Ll^^k 

k 

+ X! i^Q kPP + LR,kPR)L]^,k 

k 

- 2 X/ i-^y,k^s,k + LR^kLQ^k)pp ^ 

k 

Ppi^s,k^P.k + L^Q^f^Lp^u) 

k 

- 2Y^^^P^kLp.k + L^kLp^k, pr] ■ 

k 

Since T-Ls is invariant, the evolution of the state's i?-block turns 
out to be decoupled from the rest: in fact, by substituting the 
invariance conditions Q into ( [T6] l, we find: 



IiRL[p\IVR = ~^[Hr,Pr] +YLR^kPRL^R,k 



I {^P,k^P>k + L^R,kLR,k:PR]- 



(17) 



Now let Pr ~ vec(n/{/9n^). By exploiting the composition 
property ( [l4] i, we have: 

PR = ^RPR, 

where Cp is the map defined in ( [TSj i. 

Suppose that T-Ls is not attractive. By Theorem [T] the 
dynamics must then admit an invariant state with support 
on T-Lr. This implies that Cp has at least one non-trivial 
steady state, corresponding to a zero eigenvalue. To prove 
the converse, suppose that (0, vec(X)) is an eigenpair of Cp. 
Clearly, X ^ by definition of eigenvector. Then, any initial 
state p such that its i?-block, pp, has non- vanishing 

projection on X cannot converge to Jsi'H), and thus Hs is 
not attractive. Since D{H) contains a set of generators for 



B{H) (e.g. the pure states), there is at least one state such 
that trace(pX) ^ 0. ■ 

Building on Theorem |3] the following Corollary gives a 
bound on the asymptotic convergence speed to an attractive 
subspace, based on the modal analysis of LTI systems: 

Corollary 3 (Asymptotic convergence speed): Consider a 
QDS on n = ns (S Up, and let Us be a GAS subspace 
for the given QDS generator Then any state p e ©("H) 
converges asymptotically to a state with support only on T-Ls 
at least as fast as fce'^"*, where is a constant depending on 
the initial condition and Aq is given by: 



Ao = max{9^e(A) | A e M^r))- 



(18) 



Note that, in the case of one-dimensional Hs, the "slowest" 
eigenvalue Aq is also the smallest Lyapunov exponent of the 
dynamical system ([T]i. 

B. Connected basins approach 



Recall that the DID derived in Section III-A is a decompo- 
sition of the systems's Hilbert space in orthogonal subspaces: 



n^ns®n 



(1) 



(2) 



(9) 



By looking at the block structure of the matrices H, Lk 
induced by the DID, we can classify each basin depending on 
how it is dynamically connected to the preceding one in the 
DID. Beside Hs, which is assumed to be globally attractive 
and we term the collector basin, let us consider a basin 



'Hp = o'"' generally. Hp ~ Ti 



We can distinguish the following three possibilities for T-Lp ■ 
A. Transition basin: This allows a one-way connection 
from H^j) to y.^ ^\ when the following conditions are 
satisfied: 

=OVfc, 



ipj,^'' 7^ for some k, 



in addition to the invariance condition 



iH 



2 ^ 

k 



r(i-l)t r(«-l) 
^SM ^P,k 



(19) 



In other words, Lp^, enacts a directed probability flow 
towards the beginning of the DID: states with support 
on H^rp are "pulled" towards 

Mixing basin: This allows for the dynamical connection 
between the subspaces to be bi-directional, which occurs 
in the following cases, or types: 

I. As in the transition basin, but with 



9 E A 



'P,k 



2. In the generic case, when both L 



P,k 



+ 0, 



-'Q,k' 



^ for some fc, k'\ 



When L^fc^' = OVfc, L^' 



,fe 



^ 0, for some k. 

Circulation basin: In this case, Lp^ = and Lq ^ = 
for all fc, and thus the dynamical connection is 
enacted solely by the Hamiltonian block Hp. Not only 
is the dynamical connection bi-directional, but it is also 
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"symmetric" and, in the absence of dynamics internal to 
the subspaces 'H^\'H^^^\ and further connections, the 
state would keep "circulating" between the subspaces. 
How is this related to the speed of convergence? Let us 



consider a pair of basins "Hy 



(i) 

Hy , and let us try to 

(i) 



investigate how fast a state with support only in Hj, can 
"flow" towards in a worst case scenario. The answer 

depends on the dynamical connections, that is, the kind of 
basin the state is in, and the required speed can in each case 
be estimated as follows: 

(i) Transition basin, type-1 and type-2 mixing basins: The 
exact speed of "exit" from a transition basin has 
been calculated in |5 |, and reads 



Up) 



trace 



L 



P.k 



(20) 



which in the worst case scenario corresponds to the 

-1) 

Jk ^P,k ^P,k ' 



minimum eigenvalue of Lp ,}^^ 



.in{A|AEsp(5:L(^-)t^(^-))} 



The same quantity works as an estimate for the mixing 
basin of type-1 and-2, since to first order, with a state 
with support on T-6j) alone, the invariance condition ([19]) 



and the Lq ^ blocks have no influence. Let us recall 
that in order to have a GAS subspace, at least T-6p has 
to be a transition basin, 
(ii) Mixing basin of type-3 and circulation basin: When 



L 



(i-i) 

P,k 



for all k, and Lq j}^ ^ 0, for some k, the 



(i) 

exit from is determined by the Hamiltonian, and 
hence it is a second order effect. Therefore, it is not 
possible to estimate the "transfer speed" based on the 
derivative of the trace as we did above, since the latter 
would be zero. Similarly, when the dynamical connection 
is enacted only by the Hamiltonian block Hp, then again, 
it is a second order effect. 

Let us restrict our attention to the relevant subspace. 



n 



n 



(i-i) 



n 



T ' 



and write the restricted 



Hamiltonian in block-form: 



n 









l.i) 


rr(*-l)t 
P 





the DID and it is such that 



We can always find a unitary change of basis Urp 

(i) 

Uj. that preserves 
diag(si 



s^-^) > 0, with E^;-^' 



gular values of iJp 



being the diagonal matrix of the sin- 
in decreasing order Then the 
effect of the off-diagonal blocks is to couple pairs of 
the new basis vectors in ^\'H^^ generating simple 
rotations of the form e^''*j'^^*. Hence, any state in H^j) 



(i — l) 

will rotate towards Jij, as a (generally time-varying, 
due to the diagonal blocks of the Hamiltonian) com- 
bination of cosines, J^k^kit) cos{skt), for appropriate 
coefficients. The required estimate in this case can be 
obtained by comparing the speed of transfer induced by 



the Hamiltonian coupUng to the exponential decay in the 
previous case. When the noise action is dominant, can 
be thought as 1/T, with T being the "decoherence time 
scale" needed for the associated decaying exponential to 
reach the value e . Comparing with the action of Hp 
by letting s,; — mmj {sj} then yields: 

7/^ = Si/ arccos(e^^). 

This formula does not takes into account the effect of the 
diagonal block of the Hamiltonian, whose influence and 



optimization will be studied in Subsection IV-C 

It is worth remarking that the "transfer" is monotone in 
case (i), whereas in case (ii) it is so only in an initial time 
interval. Once we establish an estimate for all the transition 
speeds, we can think of the slowest speed, call it 7min, as 
the "bottleneck" in order to attain fast convergence. If, in 
particular, 7,nin = for a certain i, since the latter is not 
affected by the Hamiltonian, it provides a fundamental limit 
to the attainable convergence speed given purely Hamiltonian 
control resources. In fact, imagine the system to be prepared 
in a state of , the last basin in the DID, and follow its flow 
towards the attractive collector basin Hs- In the worst case, it 
will have to move through the slowest possible connection, 
associated to 7min- Conversely, connections enacted by the 
Hamiltonian can in principle be optimized, following the 
design prescriptions we shall outline below. 

In spite of giving only qualitative indications, the advantage 
of this approach is twofold: (i) Estimating the transition speed 
between basins is, in most practical situations, both simpler 
and more efficient than deriving closed-form expressions for 
the eigenvalues of the generator; (ii) Unlike the system- 
theoretic approach, it gives an immediate idea of which control 
parameters have a role in influencing the speed of convergence. 



C. Tuning the convergence speed via Hamiltonian control 

It is well known that the interplay between dissipative 
and Hamiltonian dynamics is critical for controllability [3], 
invariance, asymptotic stability and noiselessness |4|, |31, as 
well as for purity dynamics [101 . By recalling the definition of 
Cfi given in Eq. ( [l3] i. Corollary [3] highlights that not only can 
the Hamiltonian have a key role in determining the stability 
of a given state, but it can also influence significantly the 
convergence speed. In order to gain concrete insight, we first 
study a simple prototypical example. 

Example 2: Consider a three-dimensional system driven by 
a generator of the form ([T]i, with operators H, L that with 
respect to the (unique, in this case) DID basis {|s), |ri), |r2)} 
have the following form: 





- J. 





" 




' 


£ 


" 


H = 





A 


n 


, L = 
















n 

















(21) 



It is easy to show, by recalling Proposition 1, that pd — \s){s\ 
is invariant, and that any choice of il,£ 7^ also renders pd 
GAS. We can study the eigenvalues of £u as defined in ([T5]l 
and invoke Theorem [3] 




Fig. 1. Convergence speed to the GAS target state = \s){s\ for the 
3-level QDS in Example 2 as a function of the time-independent Hamiltonian 
control parameters A and fi. 



Without loss of generality, let us set £ — 1 so that L ~ 
\s) (ri|, and assume that A, ft are positive real numbers, which 
makes all the relevant matrices to be real. Let Ho := ipi^. 



We can then rewrite 



£b. = ® Ir + h: 



(22) 



where = ±iHj^ — no/2. Let be the eigenvalues of 
R'^,R~. Given the tensor structure of Cr, the eigenvalues of 
( |22] l are simply aij = X'^ + XJ , with i,j — 1,2. The real 
parts of the can be explicitly computed: 

^Re(a,,) = [-1/2,-1/2, -1/2 ± l/2mc{VT)], 

where L = 1 - + i A - 4f]2. The behavior of Aq = -1/2 + 
l/2fHe(A/r) is depicted in Figure [l] Two features are apparent: 
Higher values of lead to faster convergence, whereas higher 
values of A slow down convergence. The optimal scenario 
(|Ao| = 1/2) is attained for A = 0. 

The above observations are instances of a general behavior 
of the asymptotic convergence speed, when the Hamiltonian 
provides a "critical" dynamical connection between two sub- 
spaces. Specifically, the off-diagonal part of Hj^ in ( |2T| is 
necessary to make pd GAS, connecting the basins associated 
to |ri) and |r2). Nonetheless, the diagonal elements of H 
also have a key role: their value influences the positioning of 
the energy eigenvectors, which by definition are the system 
states that are not affected by the Hamiltonian action. Intu- 
itively, under the action of H alone, all other states "precess" 
unitarily around the energy eigenvectors, hence the closer the 
eigenvectors of H are to the the states we aim to destabilize, 
the weaker the destabilizing action will be. 

A way to make this picture more precise is to recall that 
Pd — \s){s\ is invariant, and that the basin associated to |ri) is 
directly connected to pd by dissipation. Thus, in order to make 
Pd GAS we only need to destabilize |r2) using H. Consider 
the action of H restricted to Hr = span(|ri), |7'2)). 

The Hamiltonian's i?-block spectrum is given by 



MHr) = { 



A±\/A2 



Fig. 2. Energy level configuration of the 4-level optical system discussed 
in Section V.A. Three degenerate stable states are coupled to an excited state 
trough separate laser fields with a common detuning A and amplitude Qi. 



with the correspondent eigenvectors: 



l±) 



2n 



1 



(24) 



2n 

Decreasing A, the eigenvectors of Hr tend to (|ri) ± 
|r2))/\/2, and if A = 0, the state |r2), which is unaffected 
by the noise action, is rotated right into \ri) after half-cycle. 
Physically, this behavior simply follows from mapping the 
dynamics within the i?-block to a Rabi problem (in the ap- 
propriate rotating frame), the condition A = corresponding 
to resonant driving ifTTl . 

Beyond the specific example, our analysis suggests two 
guiding principles for enhancing the speed of convergence via 
(time-independent) Hamiltonian design. Specifically, one can: 

• Augment the dynamical connection induced by the 
Hamiltonian by larger off-diagonal couplings; 

• Position the eigenvectors of the Hamiltonian as close as 
possible to balanced superposition of the state(s) to be 
destabilized and the target one(s). 

V. Applications 

In this Section, we analyze three examples that are directly 
inspired by physical applications, with the goal of demonstrat- 
ing how the control-theoretic tools and principles developed 
thus far are useful to tackle stabilization problems in realistic 
quantum-engineering settings. 

A. Attractive decoherence-free subspace in an optical system 

Consider first the quantum-optical setting investigated in 
irTSl, where Lyapunov control is exploited in order to drive 
a dissipative four-level system into a decoherence free sub- 
space (DFS). A schematic representation of the relevant QDS 
dynamics is depicted in Figure |2] Three (degenerate) stable 
ground states, |i)i=i,2.3, are coupled to an unstable excited 
state |e) through three separate laser fields characterized by 
the coupling constants fi^, i = 1,2,3. In a frame that rotates 
with the (common) laser frequency, the Hamiltonian reads 



H = A|e)(e|+^a(|e)(z| + |z)(e|), 



(25) 



(23) 



where A denotes the detuning from resonance. The decay of 
the excited state to the stable states is a Markovian process 
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characterized by decay rates 7^, i = 1,2,3. The relevant 
Lindblad operators are thus given by the atomic lowering 
operators Li — * = 1, 2, 3. 

Let the coupling factors il, be parameterized as 

r^i = O sin 9 cos <j), 

^2 = ^ sin 9 sin (f>, (26) 
= i7 cos 9, 

where > and 0<9<Tr, 0<(j)< 2-k, respectively. It is 
known (see e.g. |18| and references therein) that the resulting 
generator admits a DFS "Hdfs spanned by the following 
orthonormal basis: 



\di) = -sin0|l) + cos0|2), 
1^2) = cos6'(cos(/)|l) + sin(/)|2)) 



sin6l|3). 



(27) 



provided that 9 ^ kn and 7^^ + fcvr. In order to formally 
establish that this DFS is also GAS for almost all choices of 
the QDS parameters A,57i,7i, we construct the DID starting 



from Hs = 'Hdfs, and obtaining T-6p — span{|e)}, Ti}^ 
% (Hdfs ffi T-L^rp). The corresponding matrix representation 
of the Hamiltonian and noise operators becomes: 



,(2) 



H 



L. = 

















































A 


n' 
















n' 























sin (f) 


" 












cos 


(f) cos 9 





















•1 










VTiIcos 


(j) sin 6*1 


_ 
















1/72 cos 



















2 cos 9 sin (j) 

























^/72sign 


(cos (/)) sin (f>\ 


sin 61 1 






















' 












\/73 sin 
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V-sign 


(cos (/) sin 9) cos 9 






where sign(x) is the sign function and fi' := 
nsign(sin0cos(/)). By Proposition [T] it follows that 
■Hdfs is invariant. Furthermore, the vectorized map governing 
the evolution of the state's i?-block in ( [T5] l has the form: 



1^1 



in' 



-in' 



in' 



iA 



in' 




iA 



'in' 



■ 

in' 

-~in' 





Then, by Theorem [3] a sufficient and necessary condition for 
^DFS to be GAS is that the characteristic polynomial of Cr, 
A^ (s), has no zero root. Explicit computation yields: 



,7. 



(28) 



which clearly vanishes in the trivial cases where = \/i 
or = 0. Furthermore, there exist only isolated points in 
the parameter space such that vanishes, namely those 




20 



Q 



Fig. 3. Asymptotic convergence speed to the target DFS as a function of 
the parameters A and f2. We fixed 7^ = 0.9 and 8 = n/A and (j> = 3/4tt. 
The value of Aq is computed by means of Eq. jl8) . 



with only one 7^ 7^ 0, and the corresponding i7, — n (recall 
(p6|)). For all the other choices of the parameters, Hdfs is 
attractive by Theorem |3] Notice that the Hamiltonian off- 
diagonal elements are strictly necessary for this DFS to be 
attractive, whereas the detuning parameter does not play a 
role in determining stability. As we anticipated in the previous 
section, however, the latter may significantly influence the 
convergence speed to the DFS for a relevant set in the 
parameters-space. 

In Figure [3] we graph the value of Ao given in Eq. ( fTSj ) 
as a function of A and n, for fixed representative values of 
7i, 9, and (p. As in Example 2, small coupling n as well as 
high detuning A slow down the convergence, independently 
of -fi. That a strong coupling yields faster convergence directly 
reflects the fact that this is fundamental to break the invariance 
of the subspace span(|r2)). In order to elucidate the effect 
of the detuning, consider again the spectrum of Hn, which 
is given by Eqs. ( p3] l- ( |24] i. In the limit of 57 ^ 0, there 
exists an eigenvalue A — > 0, and the same holds for A — > cxi. 
Furthermore, the corresponding eigenvector tends to \r2) in 
each of these two limits. Thus, increasing the detuning can 
mimic a decrease in the coupling strength, and vice-versa. 

Notice that, unlike in Example 1, there is in this case a non- 
trivial dissipative effect linking |e) to \r2), represented by the 
non-zero i?-blocks of the Li's, however our design principles 
still apply. In fact, the Hamiltonian's off-diagonal terms are 
still necessary for "Hdfs to be GAS. 



B. Dissipative entanglement generation 

The system analyzed in Example 1 is a special instance 
of a recently proposed scheme |9| for generating (nearly) 
maximal entanglement between two identical non-interacting 
atoms by exploiting the interplay between collective decay 
and Hamiltonian tuning. Assume that the two atoms are 
trapped in a strongly damped cavity and the detuning of the 
atomic transition frequencies oji, i — 1,2, from the cavity 
field frequency a; can be arranged to be symmetric, that is, 
wi — = A = —{uj2 — Under appropriate assumptions 
|9|, the atomic dynamics is then governed by a QDS of the 
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form (|9|, where Eqs. ([TOli-([TT]l are generalized as follows: 

L ^ ^ (a+ (E) I + I ^ (7+), 

H= y{<^z I - I <8) cr^) + a{a^ (g) I + I (g)(T^), 

and where, without loss of generality, the parameters 7, A, a 
may be taken to be non-negative. The QDS still admits an 
invariant pure state pd = \tjjo){tfjo\, which now depends on the 
Hamiltonian parameters A, a: 

|^o> = ^(A|00)-a(|01)-|10))), n = VA^ + 2a^. 

The DID construction works as in Example 1 (where 7 = A = 
a = 1), except for the fact that while \tp2) are defined 
in the same way as in ( [T2j l, the explicit form of fourth basis 
state ( [T3] l is modified as follows: 

1^-3) --^{- 2«|00) ~ A(|10) - |01))). 

Therefore, in matrix representation with respect to the DID 

basis {\^o), 1-4^2), IV'i), IV'3)}, we get: 





' 







' 


L = 









































r 











H = V2 




_ 




a 
n 

V2 


a 




Si 






The entangled state \tpo) is thus GAS. Given the structure of 
the above matrices, the following conclusions can additionally 
be drawn. First, the bottleneck to the convergence speed is 
determined by the element L12, more precisely by the square 
of the ratio \/2A/Q,, see (|20|. Assuming that A ^ a, the 
latter is (approximately) linear with the detuning. This has 
two implications: on the one hand, the convergence speed 
decreases (quadratically) when A — 0. On the other hand, 
a non-zero detuning is necessary for GAS to be ensured in 
the first place: for A = 0, the maximally entangled pure 
state ps = \ips){iJs\, with [V's) = ^jdOl) - |10)), cannot be 
perfectly stabilized. Likewise, although the parameter a plays 
no key role in determining GAS, a non-zero a is nevertheless 
fundamental in order for the asymptotically stable state I^/jq) 
to be entangled. 

C. State preparation in coupled electron-nuclear systems 

We consider a bipartite quantum system composed by 
nuclear and electronic degrees of freedom, which is motivated 
by the well-known Nitrogen- Vacancy (NV) defect center in 
diamond Ull, lHO), GT], While in reality both the 

electronic and nuclear spins (for ^^N isotopes) are spin- 
1 (three-dimensional) systems, we begin by discussing a 
reduced description which is common when the control field 
can address only selected transitions between two of the 
three physical levels. The full three-level system will then be 
discussed at the end of the section. 



1 ) Reduced model: Let both the nuclear and the electronic 
degrees of freedom be described as spin 1/2 particles. In 
addition, assume that the electronic state can transition from 
its energy ground state to an excited state through optical 
pumping in a spin-preserving fashion. The decay from the 
excited state, on the other hand, can be either spin preserving 
or temporarily populate a metastatic state from which the 
electronic spin decays only to the spin state of lower en- 
ergy [|23l. We describe the overall optically -pumped dynamics 
of the NV system by constructing a QDS generator A basis 
for the reduced system's state space is given by the eight states 

\Eel,Sel) (E) \sn) = \Eel,Sel,SN), 

where the first tensor factor describes the electronic degrees 
of freedom, that are specified by the energy levels E^i — g,e, 
and the electron spin Sei = 0,1 (corresponding to the spin 
pointing up or down, respectively), while the second factor 
refers to the nuclear spin and can take the values — 0,1. 
To these states one has to add the two states belonging to 
the metastable energy level, denoted by \ms) \sn), with 
Stv as before. Notice that a "passage" through the metastable 
state erases the information on the electron spin, while it does 
preserves the nucleus spin state. 

The Hamiltonian for the coupled system is of the form 
Htot = Hg + He, where the excited-state Hamiltonian He and 
the ground-state Hamiltonian Hg share the following structure: 

+B{geiS^®lN+9nUi®S,) (29) 

■^-Y^^'-' ®S^+Sy®Sy + 2S, ® S,). 

Here, Sx.y = (''x,y, are the standard 2x2 Pauli matrices on 
the relevant subspace, while = — cr^) is a pseudo-spirj^j 
and Dg e, ^g,e, Q are fixed parameters. In particular, Ag^e will 
play a key role in our analysis, determining the strength of 
the Hamiltonian (hyperfine) interaction between the electronic 
and the nuclear degrees of freedom. On the other hand, B 
represents the intensity of the applied static magnetic field 
along the z-axis, and can be thought as the control parameter. 

In order to describe the non-Hamiltonian part of the evo- 
lution we employ a phenomenological model, using Lindblad 
terms with jump-type operators and associated pumping and 
decay rates. The relevant transitions are represented by the 
operators (including the decay rates) below: since they leave 
the nuclear degrees of freedom unaltered, they act as the 
identity operator on that tensor factor. Specifically: 

Li = 77^ |g,0)(e,0| ® Iat, 

L2 = ^7^ |g, l)(e, 1| (g) Iat, 

L-s = 77m |ms)(e, 1| (g) Ijv, .^O) 

L4 = ^70 I5, 0)(ms| g) Ijv, 

L5 - |e,0)(5,0|®lAr, 

Le = |e, (g) Iat. 

The first four operators describe the decays, with associated 
rates 7^, 7m70j whereas the last two operators correspond to 

'This different definition follows from the implemented reduction from 
a three- to a two- level system: specifically, we consider only |0, — 1) and 
neglect |1), and further map the states 0—^0 and — 1 — > 1. 
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the optical-pumping action on the electron, with a rate 7p. It 
is easy to check by inspection that the subspace 



Us :=span{|e,0,0),|5,0,0)} 



is invariant for the dissipative part of the dynamics: we next 
establish it is also GAS, and analyze the dynamical structure 
induced by the DID. 

a) Convergence analysis: Following the procedure pre- 



the available options for Hamiltonian tuning. Explicitly: 



sented in Sec. III-A we can prove that Us is attractive. This 
is of key interest in the study of NV-centers as a physical 
platform for solid-state quantum information processing. In 
fact, it corresponds to the ability to perfectly polarize the joint 
spin state of the electron-nucleus systerrj^] The proposed DID 
algorithm runs to completion (in seven iterations), with the 
following basin decomposition as output: 



(1) 



(7) 



where 



-1,(1) 

rirp 


= span{ 


ms,0)}, 


-1,(2) 
rirp 


= span{ 


e,l,0)}, 


-1,(3) 
rirp 


= span{ 


5,1,0)}, 


-1,(4) 
rirp 


= span{ 


e,0,l),|g,0,l)} 


rirp 


= span{ 


ms,l)}, 


-1,(6) 
rirp 


= span{ 


e,l,l)}, 


nj{r) 
rirp 


= span{ 


5,1,1)}- 



We do not report the block form of every operator, both for the 
sake of brevity and because it would not be very informative: 
we report instead the Dynamical Connection Matrix (DCM) 
associated to the above DID. The latter is simply defined as 
C :— Htot + Sfe -^fc, ^i'^h matrices represented in the 

energy eigenbasis, ordered consistently with the DID. While 
in general the DCM would not provide sufficient information 
to determine the invariance of the various subspaces due to the 
fine-tuning conditions derived in (|4]), having only of jump- 
type (or creation/annihilation) greatly simplifies the analysis. 
In fact, it is easy to see that a non-zero entry Cij due to some 
Lk means that the j-th state of the basis is attracted towards 
the i-th one. Thus, the DCM gives a compact representation 
of the dynamical connections between the basins, pointing to 



-Note that non-spin conserving decay and pumping processes t23J , which 
we neglected here, limit in practice the achievable polarization. 
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with he^g := De^g 
and h„ 



9eiB, h'^ := De^g - {gel + gn)B 



A, 



e,g ^e,g \yel 

e.g, tiiiva /!,„ .— Q — gnB. By definition, the block division 
(highlighted by the solid lines) is consistent with the DID, 
and all the empty blocks are zero. In some sense, the DCM 
can be seen as a "graph-based" visualization, similar to those 
commonly use to determine the controllability properties of 
closed quantum systems ll24l . ||251 . 

Since in our problem g^ <C gei, the diagonal entries in the 
Hamiltonian that are most influenced by the control parameter 
B are hg^e and h'g ^. All the other entries of the DCM are (for 
typical values of the parameters) either independent or only 
very weakly dependent on B. 

It is immediate to see that the 7p,7m,7o blocks establish 
dynamical connections between all the neighboring basins, 
with the exception of Hrp which is connected by the (B- 
independent) Hamiltonian elements Ae, Ag to Hrp , Hj^ . The 
DCM also confirms the fact that is invariant, since its first 
colu mn, ex cept the top block, is zero. In the terminology of 



Sec. 



IV-B 



"H^' is the only transition basin, is the only 



/(4) 



circulant basin, and all the other basins are mixing basins. It 
is worth remarking that any choice of the control parameter 
B ensures GAS. By inspection of the DCM, one finds that 
the bottleneck in the noise-induced connections between the 
basins is determined by the 70, 7p parameters. Since the latter 
are not affected by the control parameter, the minimum of 
those rates will determine the fundamental limit to the speed 
of convergence to Hs in our setting. 

b) Optimizing the convergence speed: The only transi- 
tions which are significantly influenced by B are the ones 
connecting 'H^'' to and By appropriately choosing 
B one can reduce the norm of he or hg to zero, mimicking 

"resonance" condition of Example 2. Assume that, as in the 

(2) 

physical system, Ae > Ag. Considering that Hj^ , associ- 
ated to he, is coupled to with the largest off -diagonal 
Hamiltonian term (Ae) and it is closer to T-ls in the DID, we 
expect that the best performance will be obtained by imposing 
he = 0, that is, by setting B = De/ gel- 

The above qualitative analysis is confirmed by numerically 
computing the exact asymptotic convergence speed, Eq. ([TSj. 
The behavior as a function of B is depicted in Figure [4] It is 
immediate to notice that the maximum speed is indeed limited 
by the lowest decay rate, that is, the lifetime 70 of the singlet 
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3.5 




Control field, B, Gauss 

Fig. 4. Asymptotic convergence speed to Hs as a function of tlie control 
parameter B for an NV-center. The blue (solid) curve is relative to the system 
with the metastable state, while the red (dashed) one is relative to a simplified 
system where the transition through the metastable state has been incorporated 
in a single decay operator L with rate 70 (see text). Typical values for NV- 
centers are: De = 1420 MHz, Dg = 2870MHz, Q = 4.945MHz, = 40 
MHz, Ag = 2.2MHz and g^i = 2.8MHz/G, g„ = 3.08 X 10-"'MHz/G. 
We used decay rates 7^ = 77MHz, MHz, -jm = 33MHz, 70 = 3.3MHz, 
and optical-pumping rate 7p = 70MHz. With these values, h'g ^ hg = 
2870 - 2.8B MHz, h'^^he = 1420 - 2.8B MHz and /i„ ^ 4.945 MHz. 



State with our choice of parameters. The maximum is attained 
for near-resonance control values, although exact resonance, 
B — De/gei, is actually not required. The second (lower) 
maximum correspond to the weaker resonance that is attained 
by choosing B so that hg = 0. Physically, tuning the control 
parameter such that /le = corresponds, in our reduced model, 
the excited-state "level anti-crossing" (LAC) condition that has 
been experimentally demonstrated in |fT9l . 

In Figure [4j we also plot (dashed line) the speed of conver- 
gence of a simplified system where the transition through the 
metastable state and its decay to the ground state have been 
incorporated in a single decay operator L = L4L3, with a 
rate equal to 70. This may seem convenient, since once the 
decay to the metastable level has occurred, the only possible 
evolution is a further decay into |g, 1, 1). However, by doing 
this the overall convergence speed appears to be substantially 
diminished, although still in qualitative agreement with the 
predicted behavior (presence of the two maxima, and speed 
lower than the minimum decay rate). The reason lies in the 
fact that in this reduced model, the T-6p , "H^"* transition basins 
become (part of) mixing basins, thus the non-polarizing decay 
and the Hamiltonian can directly influence (slowing down) the 
decay dynamics associated to L. 

2) Extended model and practical stabilization: A phys- 
ically more accurate description of the NV-center system 
requires representing both the electron and nucleus subsystems 
as three-level, spin-1 systems. In this case, a basis for full state 



space is given by the 21 states 

\Eei, Sei) (E) \sn), \ms) (g> \sn), 

where now the electron spin can take values Sei = 1,0,-1 
and, similarly, the nucleus spin is labeled by sn = 1,0,-1. 
The Hamiltonian is of the form: 

Hg^e = Dg^eS^_ (g) In + Q lei (8) si 

+B{gelS^(E)lN+gnlel<S)S^) (31) 
+ Ag,^{S^ (E)S^ + Sy (E)Sy+S,(E) S,), 

where now Sx.y.z are the angular momentum operators for 
spin-1. The non-Hamiltonian part of the evolution is given by 
the operators in Eq. ( pO] ) (where now 0, 1 correspond to the 
spin-1 eigenstates), to which one needs to add the dissipation 
channels associated to Lindblad operators: 

L7 = ^\g,-l){e,-l\<»lN, 
Ls = V^\ms){e,-l\(E)lN, 
Lg = y%\e,-l){g,~l\(S)lN- 

Similar to the spin 1/2 case, the dissipative dynamics alone 
would render GAS the subspace associated to electronic spin 
Sei = 0, that is Heifi = span{|e, 0, s„), \g, 0, s„)}. Thus, one 
may hope that Hs = span{|e, 0, 0), |g, 0, 0)} could still be 
GAS under the full dynamics. We avoid reporting the whole 
DID and the DCM structure, since that would be cumbersome 
and unnecessary to our scope: the main conclusion is that 
in this case nuclear spin polarization cannot be perfectly 
attained. While the hyperfine interaction components of the 
Hamiltonian still effectively connect the subspaces with nu- 
clear spin 0, 1, they also have a detrimental effect: T-Ls is no 
longer invariant. In fact, represented in the DID basis, the 
Hamiltonian has the following form: 
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The presence of Ag, A^ in the Hp block suffices to destabilize 
lis, by causing the invariance conditions in Q to be violated. 
However, these terms are relatively small compared to the 
dominant ones, allowing for a practical stabilization attempt. 
Following the approach outlined in Sec. |III-C| we neglect the 
Hp term and proceed with the analysis and the convergence- 
speed tuning. Again, the optimal speed condition is attained for 
S in a nearly-resonant LAC condition. By means of numerical 
computation, one can then show that the system still admits a 
unique, and hence attractive, steady state (which in this case 
is mixed) and that the latter is close to the desired subspace. 
In fact, with the same parameters we employed in the spin- 
1/2 example, one can ensure asymptotic preparation of a state 
with polarized ~ 1 spin with a fidelity of about 97%. 
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VI. Conclusions 

We have developed a framework for analyzing global 
asymptotic stabilization of a target pure state or subspace 
(including practical stabilization when exact stabilization can- 
not be attained) for finite-dimensional Markovian semigroup 
driven by time-independent Hamiltonian controls. A key tool 
for verifying stability properties under a given semigroup 
generator is provided by a canonical state-space decomposi- 
tion into orthogonal subspaces, uniquely determined by the 
effective Hamiltonian and Lindblad operators (the DID), for 
which we have provided a constructive algorithm and an 
enhanced version that can accommodate control constraints. 
In the second part of the work, we have tackled the important 
practical problem of characterizing the speed of convergence 
to the target stable manifold and the extent to which it can be 
manipulated by Hamiltonian control parameters. A quantitative 
system-theoretic lower bound on the attainable speed has 
been complemented by a hydraulic-inspired connected-basin 
approach which builds directly on the DID and, while quali- 
tative, offers more transparent insight on the dynamical effect 
of different control knobs. In particular, such an approach 
makes it clear that even control parameters that have no effect 
on invariance and/or attractivity properties may significantly 
impact the overall convergence speed. 

While our results are applicable to a wide class of controlled 
Markovian quantum systems, a number of open problems 
and extensions remain for future investigation. From a theory 
standpoint, it could be interesting to study the speed of 
convergence to the system's equilibria or limit sets without 
a priori knowing their structure. To this end, one could 
envision to first determine the minimal collector basin of the 
dynamics, and then analyze stability with that as a target (note 
that the fact the minimal collector basin contains the target 
subspace does not ensure invariance of the latter). Likewise, 
for practical applications, an important question is whether 
similar analysis tools and design principles may be developed 
for more general classes of controls than addressed here. In 
this context, the case where a set of tunable Lindblad oper- 
ators may be applied open-loop, alone and/or in conjunction 
with time-independent Hamiltonian control, may be especially 
interesting, and potentially relevant to settings that incorporate 
engineered dissipation and dissipative gadgets, such as nuclear 
magnetic resonance ll26l or trapped-ion and optical-lattice 
quantum simulators ll27l . Il28l . 
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